home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_7 / a7_3.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  2.3 KB  |  103 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 7.3 (Recursive Trapezoidal Rule).
  9. % Section    7.3,    Recursive Rules and Romberg Integration, Page 378
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program implements the recursive trapezoidal rule.
  13.  
  14. % The function f(x) is integrated over the interval [a,b].
  15.  
  16. %    Define and store the function f(x)    in the M-file  f.m 
  17. % function y = f(x)
  18. % y = 1 + exp(-x).*sin(4*x);
  19.  
  20. delete f.m
  21. diary f.m; disp('function y = f(x)');...
  22.            disp('y = 1 + exp(-x).*sin(4*x);');...
  23. diary off;
  24.  
  25. f(0); % Remark. f.m rctrap.m grpoly.m are used for Algorithm 7.3
  26.  
  27. pause    % Press any key to continue.
  28.  
  29. clc;, clg;
  30. %    Plot f(x) over the interval [a,b].
  31.  
  32. clc; clg;
  33. a = 0;
  34. b = 1;
  35. h = (b-a)/200;
  36. x0 = a:h:b;
  37. y0 = f(x0);
  38. c = 0;
  39. d = 2;
  40. axis([a b c d]);...
  41. plot(x0,y0,'g');...
  42. hold on;...
  43. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  44. xlabel('x');...
  45. ylabel('y');...
  46. title('The curve y = f(x) over [a,b].');...
  47. grid;...
  48. axis;...
  49. hold off;...
  50. shg; pause    % Press any key to continue.
  51.  
  52. clc;
  53.  
  54. % The recursive trapezoidal rule is applied over [a,b].
  55.  
  56. %    Place the endpoints of [a,b] in  a  and  b.
  57.  
  58. %    Place the number of times to bisect [a,b] in  N.
  59.  
  60. a  = 0;
  61. b  = 1;
  62. n  = 6;
  63.  
  64. T = rctrap('f',a,b,n);
  65.  
  66. pause    % Press any key to continue.
  67.  
  68. clc;
  69. for i=0:n, J(i+1) = 2^i; end
  70. values = [J;T];
  71. Mx1 = 'The recursive trapezoidal rule results: ';
  72. Mx2 = '   # subintervals     trapezoidal rule';
  73. clc,echo off,diary output,...
  74. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(values'),diary off,echo on
  75.  
  76. pause    % Press any key to see a nice graph y = f(x).
  77.  
  78. clg;
  79. c = 0;
  80. d = 2;
  81. axis([a b c d]);...
  82. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  83. hold on;...
  84. for i=1:(n+1),
  85.   if J(i) <= 16,
  86.     grpoly(a,b,J(i),1);
  87.     hold on;...
  88.   end
  89. end;...
  90. xlabel('x');...
  91. ylabel('y');...
  92. title('The recursive trapezoidal rule.');...
  93. grid;...
  94. axis;...
  95. hold off;...
  96. shg; pause    % Press any key to continue.
  97.  
  98. Mx1 = 'The recursive trapezoidal rule was applied with ';
  99. Mx2 = [Mx1,num2str(n),' recursions.'];
  100. Mx3 = 'The approximate value of the integral is:';
  101. clc,echo off,diary output,disp(''),disp(Mx2),...
  102. disp(''),disp(Mx3),disp(T(n+1)),diary off, echo on
  103.